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Abstract 

For a rigorous quantum simulation of nonadiabatic dynamics of electrons and nuclei, knowledge of 
not only first-order but also second-order nonadiabatic couplings (NAC), is required. Here we pro- 
pose a method to efficiently calculate second-order NAC from time-dependent density functional 
theory (TDDFT), on the basis of the Casida ansatz adapted for the computation of first-order 
NAC, which has been justified in our previous work and can be shown to be valid for calculating 
second-order NAC between ground state and singly excited states within the Tamm-Dancoff ap- 
proximation. Test calculations of second-order NAC in the immediate vicinity of Jahn-Teller and 
Renner- Teller intersections show that calculation results from TDDFT, combined with modified 
linear response theory, agree well with the prediction from the Jahn-Teller / Renner- Teller models. 
Contrary to the diverging behavior of first-order NAC near all types of intersection points, the 
Cartesian components of second-order NAC are shown to be negligibly small near Renner- Teller 
glancing intersections, while they are significantly large near the Jahn-Teller conical intersections. 
Nevertheless, the components of second-order NAC can cancel each other to a large extent in Jahn- 
Teller systems, indicating the background of neglecting second-order NAC in practical dynamics 
simulations. On the other hand, it is shown that such a cancellation becomes less effective in an 
elliptic Jahn-Teller system and thus the role of second-order NAC needs to be evaluated in the 
rigorous framework. Our study shows that TDDFT is promising to provide accurate data of NAC 
for full quantum mechanical simulation of nonadiabatic processes. 

PACS numbers: 
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I. INTRODUCTION 



Nonadiabatic transitions, i.e., transitions between adiabatic states, are ubiquitous in 
physical, chemical and biological systems.-"- In recent years there has been growing interest 
in quantum mechanical study of nonadiabatic transitions,-"- which has been regarded as a 
challenging field for theorists: Although most ab initio theories are built upon the Born- 
Oppenheimer approximation to separate the nuclear and electronic degrees of freedom, this 
approximation will break down in the region where nonadiabatic transitions occur. In order 
to describe nonadiabatic processes, it is necessary to go beyond the Born-Oppenheimer ap- 
proximation and take account of nonadiabatic couplings (NAC), which is the driving force 
for nonadiabatic transition to different potential energy surfaces (PES).- Since NAC (pref- 
erentially called as first and second derivative couplings in quantum chemistry) are defined 
as matrix elements of the first and second derivatives with respect to nuclear coordinates 
between adiabatic states (many-body wavefunctions), nonadiabatic dynamics simulation 
has long been relying on wavefunction-based methods to provide the NAC data. For more 
efficient calculation of NAC, density functional methods,— especially those based on time- 
dependent density functional theory (TDDFT), have been developed in the last decade. The 
study was initiated by Chernyak and Mukamelli who proposed to perturb the ground state 
using the nuclear derivative of Hamiltonian and to compute NAC from the density response. 
This scheme was first implemented by Baer-i?- to study H3 using a real-time approach and by 
Hu et alM^ to systematically study small molecules using the frequency-space formalism 
of Casida.— To avoid the pseudopotential problem in the calculation of NAC, all-electron 
TDDFT schemes have been independently developed by Hu et al}^ and Send et alM Alter- 
natively, formulations of NAC from TDDFT has also been achieved by Tavernelli et alM^— 
using the Casida ansatz, which is promising to correctly give NACs between excited states 
within the Tamm-Dancoff approximation (TDA). A recent study by Hu et al?"^ further 
clarified relationships between different DFT/TDDFT formulations of NAC— >^'^'^ These 
NAC schemes have been applied to nonadiabatic dynamics simulations and have shown that 
TDDFT is promising for balanced cost and performance on the computation of polyatomic 
systems.— i^^"— 

So far most studies on the computation and application of NAC are focused on the first- 
order, without much discussion on the second-order. Although second-order NAC can be in 
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principle expressed by the first-order, the numerical evaluation can not be easily carried out. 
This is because not only the differentiation of first-order NAC is needed, but also a complete 
expansion in eigenstates makes the product of first-order NAC involving these states rather 
complicated. In the wavefunction-based framework, although several methods for evaluat- 
ing second-order NAC have been presented,—"— there are very few literatures on the direct 
evaluation of second-order NAC in molecular systems. Correspondingly, the practical study 
by nonadiabatic simulation seldom takes second-order NAC into consideration. A simplified 
procedure is to replace the full quantum description as the quantum-classical simulation, 
since the time evolution of the nuclear degrees of freedom is described by a Poisson bracket 
that introduces only first order derivatives.^^ On the other hand, even the full nonadiabatic 
operators, both first- and second-order NAC, are taken into consideration in the formula- 
tion, such as ab initio multiple spawning, the second-order NAC are just ignored in the 
practice .-i^>^ It is noted that second-order NAC are often found to be small by experience,- 
however, they are not the second-order item in the Taylor expansion but originated from 
the presence of the scalar Laplacian. Therefore, in contrast to the vector form of first-order 
NAC, the second order are scalars. In order to verify the validity of neglecting second-order 
NAC in nonadiabatic simulations, it is crucial to examine the behavior of second-order NAC 
when the intersection points are approached. If the similar diverging behavior as first-order 
NAC is observed, the neglect of second-order NAC needs to be critically reconsidered. 

The aim of the present study is to develop an efficient TDDFT method for the calcu- 
lation of second-order NAC, which is desired to have the same-level computational cost as 
the first-order, and then to examine the behavior of second-order NAC near intersection 
points. For the efficiency, the explicit expansion into first-order NAC should be avoided. 
We will show that this can be achieved by using the Casida ansatz adapted for the first-order 
NAC,.^ while there is no need to explicitly construct auxiliary excited-state wavefunctions.— 
Justification of our procedure can be shown within the TDA. To check if the second-order 
NAC diverge at intersection points, we carry out TDDFT calculations within modified lin- 
ear response theory^^'^- in the immediate vicinity of Jahn- Teller,— Renner- Teller— >^ and 
elliptic Jahn- Teller intersections,^^ and compare results with model analysis. It is verified 
that our TDDFT results are in good agreement with the predictions from the models. In 
the vicinity of different types of intersections different behaviors of second-order NAC are 
revealed, either in the Cartesian components (x, y and z) or as a whole scalar: The compo- 
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nents are shown to be negligibly small near Renner- Teller glancing intersections, while they 
are significantly large near the Jahn- Teller conical intersections. Nevertheless, the compo- 
nents of second-order NAC can cancel each other to a large extent in Jahn- Teller systems, 
indicating the background of neglecting second-order NAC in practical dynamics simula- 
tions. On the other hand, it is also shown that such a cancellation becomes less effective in 
an elliptic Jahn- Teller system and thus the role of second-order NAC needs to be evaluated 
in the rigorous framework. 

The present paper is organized as follows. In Sec. II, we present the formulation of 
second-order NAC from TDDFT and its extension within modified linear response theory. 
In Sec. Ill, implementation in the planewave pseudopotential framework and computational 
details are given. In Sec. IV, practical calculations on various molecular systems possessing 
Jahn- Teller, Renner- Teller or elliptic Jahn- Teller intersections are performed, and compared 
with ideal values predicted by Jahn- Teller / Renner- Teller models. In Sec. V, we conclude 
our work. 



II. FORMULATION 



A. Second-order NAC from the adapted Casida ansartz 

In the previous work of Hu et al. rigorous TDDFT formulations of first-order NAC 
have been achieved, using the Kohn-Sham matrix elements of either /i-operator— "^i^ or 
(i-operator,— i.e., 

hn = — , dn = S , (1) 

where H is the many-body Hamiltonian and R^ is the nuclear coordinate with yU representing 
X, and z components and atom index. 

The /i-matrix formulation gives first-order NAC as 

Pj f) fr 

(v&ol ^ = ui' (v&ol = uf/\lS-'/'F,. (2) 

where \l/o (^P/) is the many-body electronic wavefunction of the ground (J-th excited) state, 
and Uj is the excitation energy. Matrix elements of S and are given by 

^ ifkr — fir) {^It " ^kr) 
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and 

dH 

hija,^ = l^Pja) , (4) 

where ipia, ^ia, fia are, respectively, the orbital, eigenvalue, and occupation number for the 
i-th. KS state with spin a. Fj is the eigenvector of the Casida equation^^ 

OF, =a;|F,, (5) 

where 



^ija,klT ~ ^cJ,T^i,k^j,l {^W ~ ^krf' + {fia — fja) {^ja — £ia)Kijcr^klT\/ (fkr — fir) {^It — ^kr) 

(6) 

with K being the KS matrix of the Hartree and exchange-correlation (xc) kernel (A^^'^), 

K,,a,kir = j j drdv'iPia (r) iP,. (r) A^^^^^ (r, r') i,ur (r') i^w (r') . (7) 

The KS orbitals have been assumed to be real for simplicity. 

On the other hand, the (i-matrix formulation gives first-order NAG as 

(^o|^|v^.)=c.fdtsV^F„ (8) 

where 

d 

dija,,! = {^ia\ l^pja) = {lpia\ ■ (9) 

The (i-matrix formulation is derived from the original /i-matrix formulation, using the rela- 
tionship between the nuclear derivatives of many-body Hamiltonian and Kohn-Sham Hamil- 
tonian. It can avoid the problem of the pseudopotential approximation in reproducing the 
inelastic terms corresponding to the off-diagonal /i-matrix elements. 

It is interesting to note that the two TDDFT formulations of first-order NAG, Eqs. ([2]) 
and ([8]), give similar but subtly different expressions for the connection between TDDFT 
quantities and many-body theory, i.e., 

hlS-'^'F, = uy'{^o\K\'^i) (10) 

and 

dlS'/'F, = ui'/'{<i^o\d,\'^i), (11) 
which can be further compared with the one for the dipole operator f^, 

rtS-V2F, = a;f (vl>o|f,|vl/,), (12) 
5 



as derived by Casida for the calculation of oscillator strength.— The expression of the 
operator, Eq. f lTT]) . shows a distinct feature as it gives different powers in S and uj. It is 
reminded that Eq. (|T2l) is the basis of the Casida ansatz, in which the auxiliary many-body 
excited-state wavefunction is constructed as 

= J'-^^^F^^.^dld^.^o, (13) 

ija 

SO that 

{^o\d^\^j) = {^o\0,\^i). (14) 

Herein aj^ and dia are respectively creation and annihilation operators, and is a Slater 
determinant of occupied KS orbitals. Details regarding the Casida ansatz and the mapping 
between TDDFT quantities and many-body theory can be found in Ref . |23|] . Nevertheless, 
in order to validate Eq. f lT4|) also for 0^ = d^, the Casida ansatz need to be adapted according 
to Eq. ffTTj) in the following way, 

^1= J2 \ — T Fij,jd]„di„^o, (15) 

where = ^o- 

With the adapted Casida ansatz in hand, we can now readily derive the second-order 
NAC, assuming the similarity between first- and second-derivative operators. Defining 

K - Ij. (16) 

we can get 

(^oI^mI^/) = (^0 K ^i)- (17) 



from the adapted Casida ansatz. Since Eq. ( |T5l) is equivalent to 

fi(T> fa 

^1= E (S^/^F,)^^.^at>.^o, (18) 

ija 

further using the connection from the Casida ansatz to the mapping between TDDFT and 
many-body theory^^, we can get 

hlS'/'Fi = u;'/'{^o\b,\^i), (19) 
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i.e., 

(*o|^|^.)=c.fbtS^^^F,. (20) 

This expression shows that we can calculate second-order NAC without explicitly construct- 
ing (auxiliary) excited wavef unctions. Moreover, it is appealing that the computational cost 
of second-order NAC by this expression is at the same-level as that of the first-order. On 
the other hand, it is noted that although the derivation of Eq. ( ITT]) is rigorous, derivation 
of Eq. (I19p is not yet. The validity of the adapted Casida ansatz for the second-order NAC 
needs to be further justified. Next we show that this can be achieved within the TDA, where 
the adapted Casida ansatz becomes equivalent to the original one. 



B. Second-order NAC within the TDA 



The justification of the second-order NAC formulation can be attempted by using the 
expansion of first-order NAC to show the validity of Eq. f|T7|) . It has been shown^iiS^ that 
for those between ground state and singly excited states, it generally holds that 



d„ 



(21) 



and for those between singly excited states, the validity of the expression 



{^i\h^\'^j) = {^i\h^\^j) 



(22) 



can be justified using the TDA, where cj^^S^/^ = 1, and the two forms of auxiliary wave- 
functions become the same, i.e., = \E'/. The second-order NAC can be expanded by the 
first-order as 



^oIt^^^/) + (^ol TTFT- I*/; 



'dR„ 



dR„ 



OR,, 



^ ^0|^m)(^m|7^*/) + 7^ (*0| 7^ 1^/) 



dR 



'dR,, 



dR,, 



dR„ 



J2 (^ol d, + ^ (v&ol d, \^,) , (23) 



Ei — E„ 



dR., 
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which rigorously holds since (\l/m|^'n) = ^mn- Similarly, if we can show this orthonormalized 
condition for the auxiliary wavefunction, i.e., (^m|^n) = ^mn, we can get 




From Eq. fl2T]) and Eq. f l22|) . the identity of Eq. fl23|) and Eq. f l2^ can be justified provided 
that \E'/ = 'I'/, since we have to reconstruct the auxiliary wavefunction from to 'I'm when 
the operator is changed from to h^. This is satisfied when the TDA is valid. In the 
meanwhile, the orthonormalized condition that 

ija klr 

= y^EE (S'^'F,)J^.^ {S'/'F j)^^^6.kSji6^r = V^F\SFj (25) 

ija klr 

also holds within the TDA since F|Fj = 6fj. Therefore, Eq. ( 123|) and Eq. (!24|) become 
identical, i.e., the validity of Eq. f lT7|) is justified within the TDA. 

Further remark is on the complete expansion in Eq. (^^. As long as we only consider 
a singly excited state, this does not pose a problem since is a single Slater determinant 
and only the contributions from other singly excited states enter the expansion. 

C. Extension within TDDFT modified linear response theory: Justification of the 
Slater transition state method 

In the calculation of first-order NAC, a particular example is the case of the Slater 
transition state method for doublet systems. Billeter and Curioni^ has used the following 
expression, 

(*olg|l*,) = (ClallC), (26) 

where the {i,j) pair is the particle-hole orbitals responsible for the I-th transition, and m 
denotes the mid-excited state (Slater transition state) in which the particle-hole orbitals 
are each filled with a half electron. They have found that this expression can give accurate 
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^ia^ja^lr^kr 



results of first-order NAC between doublet states of molecules at equilibrium geometries, 
and their approach is further validated by our TDDFT modified linear response theory2ii2^ 
and also by our calculations near intersection points.— Next we will show that the extension 
of TDDFT formulation of second-order NAC within modified linear response theory, is also 
equivalent to the Slater transition state method for doublet systems. 

Within modified linear response, the excitation energy is calculated from the response of 
the mid-excited state, while other terms in the NAC formula are calculated from that of the 
pure-state configuration.— Corresponding to the mid-excited state of a doublet system, the 
adapted Casida equation, 

0™F7 = w™F7 (27) 

with the matrix element 
gives 

< = S"-e™, (29) 

since fi^=fj^ = 0.5 in the mid-excited state of a doublet system, which renders the corre- 
sponding off-diagonal elements of f2 to be zero. On the other hand, the pure state configu- 
ration in the mid-excited state, which uses the occupation number of the ground state while 
keeping other quantities of the mid-excited state, gives 

bLpSy^Ff = b-AeT. - a-'^', (30) 

due to the fact that F[j^j is practically equivalent to 1 and other components of F/ are zero. 
Therefore, 

i^ohl^i) = {u^TY^'hlA^'F^i = = (CI^IC)' (31) 
which is just the second- derivative coupling matrix element between the particle-hole or- 
bitals. As a result, the TDDFT formulation of second-order NAC in doublet systems is just 
reduced to the Slater transition state method. 

III. IMPLEMENTATION AND COMPUTATIONAL DETAILS 

The implementation of the present TDDFT method for second-order NAC is based on 
the ABINIT code,— which is a planewave pseudopotential approach. All calculations are 
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performed within adiabatic LSDA using the Teter Pade parametrization.— The TrouUier- 
Martins pseudopotentiala^ with nonhnear core correction,— generated by Khein and Allan, 
are used for various atomic species. Only the F point {k=0) is taken into consideration 
in the k point sampling, which corresponds to the use of real wavef unctions. Convergence 
parameters, such as the supercell size, number of unoccupied orbitals, and kinetic energy 
cutoff, are examined to ensure reasonably accurate results. On the basis of the previous 
implementation of modified linear response theory in ABINIT,— its extension for calculating 
second-order NAG requires almost no additional labor, since it is only necessary to construct 
the pure-state configuration from the mid-excited state, and to apply the same calculation 
procedures as ordinary linear response theory. To check the performance of our method, 
it is desired to compare TDDFT results for general atomic geometries with those from 
wavefunction-based methods, however, there are too few literatures on this aspect and the 
direct comparison is difficult. Therefore, we concentrate on evaluating second-order NAG in 
the immediate vicinity of Jahn- Teller, Renner- Teller and elliptic Jahn- Teller intersections, 
where we can directly compare our results with predictions from corresponding models. 

Finite difference method of calculating 6-matrix elements 

The calculation of 6-matrix elements is implemented in a straightforward finite-difference 
scheme, with the consideration of aligning the phases of KS orbitals, as shown by 

|?|, , _ {iJ^aiRMA^ + AR ■ e^)sgn(e+) - 2^,^(R) + ^,^{R - AR ■ e^)sgn(e-)) 

(32) 

where is the unit vector along the /i axis, sgn(^) is the sign function, i.e., 

/ -1 < 

sgn(0 = < (33) 
[ 1 if ^ > 

and 

e+ = (^,.(R)|^,.(R + AR ■ e^)), (34) 
e- = (^,^(R)|V',.(R - AR ■ e,)). (35) 

The accuracy of the above numerical differentiation scheme is checked by using different 
AR. In the practice, we choose Ai?=0.002~0.004 bohr. 
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1 3 

FIG. 1: The geometry of the X3 system as one X atom (numbered as 2) is moved on the contour 
around the mtersection point (located at O). The nuclear configuration at the intersection point 
is an equilateral triangle with C^y symmetry, corresponding to the degeneracy of the ground state 
and the first excited state. 

IV. RESULTS AND DISCUSSIONS 

In this section, we present calculation results on various molecular systems possessing 
Jahn- Teller, Renner- Teller, or elliptic Jahn- Teller intersections, where the ground state and 
the first excited state of these molecular systems are degenerate. 

A. Jahn- Teller systems 

In Table H] we list the x, y, and z components of second-order NAC in two typical Jahn- 
Teller systems: The prototype H3 molecule^"— and an alkali-metal trimer Lis.— The 
three atoms are located in the geometry of Fig. [1] in which one atom is moved on the 
contour around the intersection point. The contour radius q is chosen as 0.02 bohr, which 
is sufficiently small so as to be comparable to the condition of the Jahn- Teller model. It is 
clearly seen that at such a small q the x and y components of second-order NAC in H3 and 
Lis are significantly large: The nonzero values are in the order of 1000 bohr^^, which are 
much larger than those of first-order NAC (which are in the order of 1/g). In the meanwhile, 
both the magnitude and relative signs of TDDFT results are in good agreement with the 
Jahn- Teller model. (The ideal values of second-order NAC from the Jahn- Teller model can 
be derived from the derivatives of the first-order NAC, as shown by Appendix |Al) On the 



11 



other hand, it is noted that the z components of second-order NAC are quite small but 
nonzero, either in H3 or Lis. This is different from the zero values of z components of first- 
order NAC in X3 systems near intersection points. The Jahn- Teller model predicts that the 
X and y components of second-order NAC only depend on contour radius q and angle ^, 
while there is an additional dependence of z component on the internuclear distance r. In 
our calculations we set th-h = 1-9729 bohr and rLi_Li = 5.0 bohr respectively, therefore, 
TDDFT calculations are expected to give different z components for H3 and Lis, and the 
results seem to give such a difference: Both the magnitude and sign agree with the ideal 
values corresponding to the above internuclear distances. However, since the z components 
are quite small we need to make sure whether they are intrinsically nonzero. For this purpose 
we have made a detailed examination of the z components of second-order NAC on the three 
atoms of H3 as a function of the contour angle ^, as shown by Fig. [2j As 6' is varied from to 
180°, the z components on all atoms, although small, show clear dependences on Q and agree 
well with the Jahn- Teller model. This means that the small z components of second-order 
NAC are intrinsically nonzero and have been accurately reproduced by TDDFT. 

Another point noteworthy in Table |T] is the sum of x, y and z components. In contrary 
to the vector form of first-order NAC, second-order NAC are scalars due to the presence of 
the scalar Laplacian, therefore, only the sum of x, y and z components are meaningful in 
the nonadiabatic dynamics simulation. Table |T] shows the sum of components in both H3 
and Li3 are small, as predicted by the Jahn- Teller model. This can provide the background 
for the the neglect of second-order NAC in practical simulations.-'^'^ 

B. Renner- Teller systems 

In Table [IT] we list the x, y, and z components of second-order NAC in several typical 
Renner- Teller systems: The BH2, CH^, NH2 and H20^ molecules,—'^ which are in the 
geometry of Fig. |3] with the contour angle 6' = 0. The contour radius q is chosen as 0.1 
bohr, which is known to be sufficiently small and can be comparable to the condition of the 
Renner- Teller model. The internuclear distances are set as th-b = 2.0 bohr, th-c = 2.0 
bohr, th-n = 1-95 bohr, and rn-o = 1-85 bohr. It is interesting to see that all components 
of second-order NAC on three atoms of all molecules are negligibly small (reminding that 
the first-order NAC in Renner- Teller systems are in the order of 1/g) and agree with the 
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TABLE I: The calculated x, y and z components of second-order NAC (in bohr~^) on three atoms 
of H3 and Li3, which are at the geometry of Fig. [TJ The contour radius q is 0.02 bohr and angle 9 
is 0. The ideal values from the Jahn- Teller model, as derived in Appendix |X] and summarized in 
Table |Vl are also listed for comparison. It is noted that the z components of second-order NAC in 
the Jahn- Teller model are dependent on atomic distances and have been derived within two sets 
of parameters: the values inside the parenthesis are derived by rLi-Li = 5.0 bohr, while the others 
outside are derived by th-h = 1-9729 bohr. 



X 



y 



H, 



atom 1 
atom 2 
atom 3 



1085.88 
0.30 
-1085.30 



-1074.36 
0.00 
1073.36 



12.75 
0.00 
-12.68 



Lia 



atom 1 
atom 2 
atom 3 



1102.08 
1.44 
-1099.12 



-1090.50 
0.00 
1086.05 



4.37 
0.00 
-4.71 



Model 



atom 1 
atom 2 
atom 3 



1082.53 
0.00 
-1082.53 



-1082.53 
0.00 
1082.53 



12.67 (5.0) 
0.0 (0.0) 
-12.67 (-5.0) 



Renner- Teller model. (The ideal values of second-order NAC from the Renner- Teller model 
can be derived from the derivatives of the first-order NAC, as shown by Appendix [Bl) As 
a matter of fact, NAC of Renner- Teller system do not dependent on the contour angle 9, 
and thus the negligibly small values of second-order NAC in demonstrated systems are not 
accidental results for a specified geometry, but indicate that they are intrinsically zero. In 
connection with the nonadiabatic dynamics simulation, it is thus verified that the sum of x, 
y and z components can be regarded as zero in Renner- Teller systems. This also provides a 
background for the the neglect of second-order NAC in practical dynamics simulations.-i^i^ 
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30 60 90 120 150 180 

9 (degree) 

FIG. 2: The z components of second-order NAC on the three atoms of H3. The labels Iz, 2z and 
3z denote the z components on atom 1, 2 and 3, while lz_JT, 2z_JT and 3z_JT denote those from 
the Jahn- Teller model, respectively. 




FIG. 3: Geometry of the XH2 or XH2 system when the X atom is moved on the contour around the 
Renner- Teller intersection point (indicated by the open square) on the collinear axis. The contour, 
with radius q and angle 0, is fixed in the xy plane, which is perpendicular to the HH axis. The 
two hydrogen atoms are set to be symmetric to the plane. 

C. Elliptic Jahn- Teller system 

In Table UTTl we list TDDFT calculation results of the x, ?/, and z components of second- 
order NAC on the three atoms of NaH2, which is known as an elliptic Jahn- Teller system.— 1^ 
The three atoms are located in the geometry of Fig. H] with the contour angle Q = 60°. Other 
parameters regarding the geometry are r = 2.18 bohr and R = 3.6127 bohr, according to the 
intersection point determined in our previous work.— For an elliptic Jahn- Teller systems, 
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TABLE II: The calculated x, y and z components of second-order NAC (in bohr^^) on three atoms 
of BH2, CHg^, NH2 and H20^, which are at the geometry of Fig. [3l The contour radius q is 0.1 bohr 
and the angle 9 is 0. The ideal values from the Renner- Teller model, as derived in Appendix IB} 
are also listed for comparison. 



y 



BHo 



atom H 
atom B 



(1) 



atom H 



(2) 



0.015 
0.016 
0.015 



0.0 
0.0 
0.0 



0.0 
0.0 
0.0 



CH+ 



atom H(i) 
atom C 
atom H(2) 



-0.012 
0.013 
-0.012 



0.0 
0.0 
0.0 



0.0 
0.0 
0.0 



atom H(i) 
atom N 
atom H(2) 



0.018 
0.018 
0.018 



0.0 
0.0 
0.0 



0.0 
0.0 
0.0 



H20+ 



atom H(-i) 
atom O 
atom H(2) 



-0.0035 
0.021 
0.0014 



0.0 
0.0 
0.0 



0.0 
0.0 
0.0 



Model 



atom 1 
atom 2 
atom 3 



0.0 
0.0 

0.0 



0.0 
0.0 

0.0 



0.0 
0.0 

0.0 



the angular NAC Aq is not just in a quantized value of 0.5, but shows a strong dependence 
on the contour angle 9. Using the similar procedures in Appendix |A] but setting the angular 
NAC, Aq, as a variable rather than a constant of 0.5, we can easily get the conclusion that 
second-order NAC would depend on i.e., the slope of angular NAC with respect to 6. 
Therefore, we set 6 as 60° at which is relatively large, as revealed by our previous work.— 
It is clearly seen that under such a condition the magnitudes of x and y components become 
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FIG. 4: The geometry of the NaH2 system as the Na atom is moved on the contour around the 
conical intersection point (indicated by the open square). The nuclear configuration at the conical 
intersection point is an isoseles triangle with C2v symmetry. 



TABLE III: Calculated values of x, y and z components of second-order NAC (in bohr"-^) on the 
three atoms of NaH2, which are in the geometry shown by Fig. H] with the contour radius q=0.1 
bohr and angle 9 = 60°. 



X 



y 



atom Na 



atom H 
atom H 



(1) 

(2) 



20.54 
-542.53 
-618.23 



-110.13 
-10.94 
51.59 



1.23 
-2.89 
-5.42 



unbalanced in comparison with the Jahn- Teller systems. Meanwhile, the z components are 
still relatively small. Therefore, the sum of x, y and z components of second-order NAC 
can not be negligibly small. Regarding the role of NAC in the nonadiabatic dynamics 
simulation, it is thus suggested to include second-order NAC in the rigorous simulation of 
general molecular systems, which might possess accidental conical intersections without any 
symmetry requirements as in Jahn- Teller systems. 



V. CONCLUSION 



We have proposed an efficient TDDFT method for calculating second-order NAC between 
ground state and singly excited states, which is based on the Casida ansatz adapted for 
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first-order NAC, while the calculation procedure can be done without the need of exphcitly 
constructing auxiliary excited-state wavef unctions. Our formulation can be justified when 
the TDA is valid. Within the modified linear response theory, the TDDFT formulation is 
reduced to the Slater transition method in doublet systems. Test calculations are carried 
out in the immediate vicinity of various types of intersection points. The results are in good 
agreement with the ideal values derived from Jahn- Teller or Renner- Teller models. Contrary 
to the diverging behavior of first-order NAC near intersections, the Cartesian components 
of second-order NAC are shown to be neghgibly small near Renner-Teller intersections, 
while they are significantly large near Jahn- Teller intersections. Nevertheless, the Cartesian 
components of second-order NAC can cancel each other to a large extent in Jahn- Teller 
systems, showing the background of neglecting second-order NAC in nonadiabatic dynamics 
simulations. On the other hand, it is revealed that such a cancellation becomes less effective 
in an elliptic Jahn- Teller system and thus the role of second-order NAC needs to be evaluated 
in the rigorous framework. Finally, it is noted that the performance of TDDFT on the 
computation of second-order NAC needs to be further validated, particularly for a general 
atomic geometry, which requires reference data and remains future work. 
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FIG. 5: Configuration of a Jalin- Teller trimer in the xy plane, in which atom 2 is regarded on 
a contour with radius q and angle 9 around the intersection point (vertex D of the equilateral 
triangle). Arrows represent NAC vectors on the three atoms. 

TABLE IV: The x, y and z components of first-order NAC on the three atoms of a Jahn- Teller 



trimer, which are in the 


geometry shown by Fig. [5j 


q is the contour radius and I 


■1 is the contour 


angle. 










X 


y 


z 


atom 1 


^ cos(120° - 6) 


-^sin(120° -0) 





atom 2 


^cose 


q 





atom 3 


-^cos(60° -0) 


^ sin(60° - 9) 






Appendix A: Second-order NAC from the Jahn- Teller model 

The Jahn- Teller model describes a class of systems in which a set of nuclear coordinates 
are coupled to a two-level system consisting of the ground state and the first excited state of 
appropriate symmetry^. Figure |5] shows an arbitrary configuration of a Jahn- Teller trimer. 
When the contour radius q is sufficiently small, the angular NAC has a quantized value of ^ 
according to the Jahn- Teller model.^^ All components of first-order NAC on the three atoms 
can thus be uniquely determined, as shown by Table IIV[ 

To derive the x component of second-order NAC on atom 2, we move atom 2 in the x 
direction with a small displacement A, as shown by Fig. |6l Since the Jahn- Teller model is 
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FIG. 6: Schematic view of the derivation process of the x component of second-order NAC on atom 
2. A small displacement A is made in the x direction for atom 2. After displacement the contour 
radius is changed from q to q' , and the contour angle is from to 6' . The arrow denotes the new 
NAC vector on atom 2. 

a two-level system, we can get 

where A^^ is the x component of first-order NAC before the displacement, as listed in 
Table HVl The new x component of first-order NAC on atom 2 after the displacement, A^xl"^ 
can be determined from the new geometry as 

<^ = ^cos^', (A2) 
q 

where 

q = cos2^ + (gsin^- A)2 (A3) 

and 

ni f qcos6\ , . ^ . 

e = arccos f J . (A4) 

By taking A — )■ in Eq. (lAip . we can get 

/X , ^ , 1 /^0.5 ^, 0.5 \ 0.5 

(^0 TT^ ^i) = lim — — cos 6* cos 6* = — sm26' 

3x2 A \ q' q J q 



(A5) 



The derivation of the y component of second-order NAC on atom 2 is similar to that 
of the X component in the above, thus the detail is not shown here. Next, to derive the z 

19 




FIG. 7: Schematic view of the derivation process of the z component of 2nd-order NAG on atom 2. 
After a small displacement A is made in the z direction for atom 2, the contour radius is changed 
from q to q' , and the NAG vector (denoted by the arrow ending on atom 2) is located in the new 
atomic plane. 

component, we move atom 2 in the z direction as shown by Fig. U\ and then we can get 

Q2 d d A^'^'P 

which uses the fact that the z component of first-order NAC on atom 2 before the displace- 
ment is zero. The new geometry after the displacement gives 

^dlSp _ ^ _ pQg (y^ pQg Q,^^ (A7) 

q' q' 

where g' is the new contour radius around the vertex of the equilateral triangle in the 
new atomic plane, a is the angle between the new NAC vector and the z axis, ai is the 
angle between the new atomic plane and the z axis, and a2 is the angle between the new 
NAC vector and the projection of the z axis in the new atomic plane. Using the geometric 
relationships shown by Fig. [TJ we can get 



q = i + r"^ — 2r^r cos 



arccos ( ^ ^ 1 - 60° 
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FIG. 8: Schematic view of the derivation process of 2nd-order NAG components on atom 1, which 
is regarded as rotating around vertex A of a new equilateral triangle with side length tq. The 
corresponding contour radius and angle is q and respectively. 

«! = arcsin (7—] , (A9) 



h 



and 



2q'r4 J \r4 
The auxiliary quantities in the above equations are calculated as 



a2 = 90° — arccos I —y ) — arccos I — j . (AlO) 



r3 = \ rl + A2, n = \ rl + 



/ii = — + g cos 6*, /i2 = yhj + A'^, 
ri = -^7-2 _|_ g2 _ 2qr cos(150° — 9), 
ra = ^/r^ + - 2qr cos(210° - 9). 

By taking A — )■ 0, we can get ^ ri, ^ r2, h2 ^ hi, q' q, and 0^2 — )■ 90° — 9. Then 
Eq. (1A6I) is reduced to 



10 5 
A^o A q 



arcsin 



h2 



cos(90° -9) = — sin 9, (All) 
q v3r/2 



where we have used the fact that r ^ q. 

To derive components of second-order NAC on atom 1 and atom 3, we need not make 
displacements but can merely use the fact that the three atoms are equivalent, i.e., not only 
atom 2 can be regarded as rotating in a contour around the intersection point, other two 
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TABLE V: The x, y and z components of second-order NAC on the three atoms of a Jahn- Teller 
trimer, which are in the geometry shown by Fig. [5j q is the contour radius and is the contour 
angle. 



atom 1 ^ cos(30° + 2(9) cos(30° + 26^) ^^-^iy^ sin(6l + 120°) 

atom 2 ^sin26' -^sm2e ^^^sinO 

atom 3 sin(120° - 26) ^ sin(120° - 26) ^TFTi ^^^(^ ~ ^^0°^ 



atoms can also be taken into such a view. In Fig. [HI where atomic geometry is the same as 
in Fig. [5l atom 1 is regarded as rotating around the intersection point A with contour radius 
q and angle 9q. Here A is the vertex of a new equilateral triangle with side length tq. The 
geometric analysis gives 



+ ^2 _ 2qr cos(210° - 6), (A12) 



and 



,q^ + rl-r'^\ „ - r cos(210° - 6*) 

= 150° + arccos ^ \ ° = 150° + arccos ^ ^ ' 



(A13) 



2gro 

Using the fact that r ^ g we can easily get ro = r and Oq = 120° + 6. Replacing r and 
6 in the expression of second-order NAC components on atom 2 with rg and 6o, we can 
immediately get the results for atom 1. 

The results for atom 3 can be derived in a way similar to that for atom 1. The final 
results of second-order NAC components on all atoms are listed in Table El 



Appendix B: Second-order NAC from the Renner- Teller model 

Figure M shows an arbitrary configuration of an XY2 Reller- Teller systems, where the X 
atom is moved with a sufficiently small displacement q from the z axis, which is the seam 
of Renner- Teller intersections. According to the Renner- Teller model^, the angular NAC 
has a quantized value of 1.0 and all components of first-order NAC on three atoms can be 
determined, as shown by Table I VII Note that all y and z components are equal to zero, 
meaning that the NAC vectors are parallel to the x axis. 
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FIG. 9: Configuration of an XY2 Renner- Teller system in the yz plane. The three atoms are 
located in a geometry slightly distorted from the linear geometry, and atom 2 is regarded on a 
contour around the z axis with contour radius q. 

TABLE VI: The x, y and z components of first-order NAC on the three atoms of an XY2 Renner- 
Teller system, which is in the geometry shown by Fig. [9l g is the contour radius, while ri (r2) is 
the distance of atom 1 (atom 3) from the intersection point. 



atom 1 
atom 2 
atom 3 



_1_£2_ 

q ri+r2 

1 
1 

1 r-i 
q ri+r2 










Because the first-order NAC vectors are perpendicular to the yz plane, small movement of 
atoms in the yz plane will not alter the direction of NAC vectors and the y and z components 
of first-order NAC are kept to be zero after the displacement. Therefore, we can immediately 
conclude that y and z components of second-order NAC on all atoms are zero. 

To derive the x component of second-order NAC on atom 2, we move atom 2 in the x 
direction, as shown by Fig. [TOT a). Since the Renner- Teller model is a two-level system, we 
can get 

f)^ f) r) /Idisp _ 4 



where 



A,, = -. (B2) 



The new geometry after the displacement gives 



Af^^ = ^cose, (B3) 
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FIG. 10: Schematic view of the derivation process of the x component of second-order NAC on (a) 
atom 2 and (b) atom 1. A smah displacement A is made in the x direction for atom 2 in (a) and 
atom 1 in (b). After displacement the contour radius is changed from qto q' . Note that in (b) the 
intersection point is changed from O to O'. 



where q' = a/q'^ + and 9 = arccos(g/g'). By taking A — t- in Eq. ( IBll) . we can get 

-A2 



92 1 
(^o|^|^i) = lim 



(g2 + A2)g 



0. 



(B4) 



A^O A 

In a similar way the x component of second-order NAC on atom 1 can be derived by 
displacing atom 1 in the x direction, shown by Fig. [TOT b). as 



where 



92 



/$o| — 1^1) = - ' 



A 



dxi dxi 
1 r-2 



A 



g ri + r2 ' 
1 r' 



q' = ^/q^ + (r2 - r'^y, 

9 = arccos(g/g'), 

ri + r2 ri + r2 

1^2 = r2- —. — 7 = ^2 ■ 



r'l + r'^ v/(ri+r2)2 + A2" 
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(B5) 

(B6) 

(B7) 

(B8) 
(B9) 
(BIO) 



Taking A — )■ in Eq. ( IBSp . we can get 



, , 92 , , , 1 gVaA^ + r^A^ 

^ axf ' ' A^oA[g2A2 + g2(ri + r2)2 + r|A2)]g(ri + r2) ^ ^ 

Finally, since the derivation of x component of second-order NAC on atom 3 is essentially 
equivalent to that of atom 1, the result for atom 3 is also 0. 

In one word, all components of second-order NAC on the three atoms of the Renner- Teller 
system are equal to zero. 
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